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Abstract 

GPS  signals  can  be  used  for  positioning  satellites  in  geostationary  (GEO)  orbits,  but  the 
performance  in  this  case  is  poor,  because  very  few  pseudorange  measurements  are  available  at 
any  given  time.  This  paper  describes  a  new  method  for  improving  geostationary  satellite 
navigation  accuracy  by  using  dynamic  Two-Way  Time  Transfer  (TWTT)  measurements.  By 
directly  measuring  the  clock  error  between  the  GPS  satellite  and  the  GPS  receiver,  TWTT 
allows  meaningful  information  to  be  gathered  when  less  than  four  GPS  satellites  are  available. 
A  simulation  was  developed  in  which  satellites  in  GEO  orbits  with  GPS  receivers  onboard 
generated  a  position  with  1)  GPS  with  a  crystal  clock,  2)  GPS  with  an  onboard  atomic  clock,  3) 
GPS  with  TWTT  to  a  ground-based  atomic  clock,  and  4)  GPS  with  TWTT  to  a  ground-based 
clock  synchronized  to  GPS  time.  Bringing  an  atomic  clock  into  the  system  (Cases  2  and  3) 
resulted  in  a  21-38%  improvement  in  the  3-D  RMS  position  accuracy  over  the  standard  GPS 
case  (Case  1).  However,  using  TWTT  with  a  clocked  slaved  to  GPS  time  resulted  in  a  60%-70% 
improvement  in  3-D  RMS  positioning  accuracy.  This  level  of  performance  was  obtained  for 
TWTT  measurement  error  standard  deviations  anywhere  between  0.3  ns  to  30  ns. 


INTRODUCTION 

Dynamic  two-way  time  transfer  (TWTT)  has  recently  been  demonstrated  [1,2],  opening  up  the  possibility 
of  using  TWTT  measurements  to  improve  GPS-based  navigation  solutions  on  moving  platforms. 
Previous  research  has  shown  a  40%  improvement  in  DGPS  positioning  accuracy  when  using  TWTT  to 
synchronize  clocks  between  a  network  of  GPS  receivers  [3].  Similar  results  can  be  found  in  [4], 

The  main  objective  of  this  research  was  to  examine  the  impact  of  adding  dynamic  TWTT  measurements 
to  GPS-based  geostationary  (GEO)  satellite  positioning.  Using  GPS  to  position  satellites  in  GEO  orbits 
can  be  a  challenging  task.  Most  of  the  L-band  RF  energy  transmitted  by  the  GPS  satellites  is  aimed  at  the 
Earth,  and  only  occasionally  is  a  GEO  satellite  within  the  main  beam  of  a  GPS  satellite,  as  shown  in 
Figure  1.  As  a  result,  very  few  pseudorange  measurements  are  typically  available  for  positioning  a  GEO 
satellite.  Figure  2  shows  an  example  of  the  number  of  available  GPS  measurements  for  a  GEO  satellite 
over  a  1  -day  period.  From  this  plot,  it  is  clear  that  having  four  satellites  (which  is  required  to  obtain  a  full 
position/clock  error  solution)  is  a  relatively  rare  occurrence  for  a  GEO  satellite.  When  fewer  than  four 
measurements  are  available,  then  the  receiver  clock  error  cannot  be  estimated,  and  receiver  clock  errors 
will  affect  the  positioning  solution. 
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Figure  1.  GPS/GEO  satellite  coverage  geometry. 


Figure  2.  Example  of  the  number  of  available  pseudorange  measurements  for  a  standard 
GPS  receiver  in  orbit  on  a  GEO  satellite. 


The  thought  behind  this  research  is  that  by  using  TWTT  measurements  to  help  constrain  or  measure  the 
GPS  receiver  clock  error,  then  the  positioning  solution  should  be  improved  for  a  GEO  satellite.  TWTT 
measurements  should  eliminate  the  need  for  a  precise  clock  on  the  satellite,  because  it  would  only  need  a 
precise  reference  clock  on  the  ground.  Essentially,  using  TWTT  with  a  highly  accurate  clock  on  the 
ground  and  a  low-quality  clock  on  the  satellite  would  be  comparable  to  putting  a  highly  accurate  clock  on 
the  satellite  itself.  Even  better  performance  should  be  possible  by  using  the  TWTT  system  to  directly 
measure  the  absolute  GPS  receiver  clock  error  on  the  GEO  satellite  by  synchronizing  the  GPS  receiver 
clock  to  a  clock  slaved  to  GPS  system  time. 

A  simulation  was  developed  in  order  to  evaluate  the  effect  of  different  ways  to  use  TWTT  measurements 
to  improve  positioning  of  GEO  satellites.  This  paper  describes  the  simulation  and  summarizes  the  key 
results  obtained  by  running  the  simulation. 


512 


39th  Annual  Precise  Time  and  Time  Interval  (PTTI)  Meeting 


BACKGROUND 

Two-Way  Time  Transfer  (TWTT)  is  a  technique  in  which  signals  are  simultaneously  exchanged  between 
two  users  to  measure  their  relative  clock  offsets.  If  the  paths  between  the  two  users  are  reciprocal,  the 
delays  cancel  and  the  difference  between  the  two  clocks  is  half  the  difference  in  time-interval  counter 
readings  [5]. 

TWTT  can  be  performed  in  both  static  and  dynamic  modes.  Static  TWTT  uses  two  or  more  transceivers 
whose  positions  on  Earth  are  held  constant  during  the  transmission  and  reception  of  the  measurement 
signals.  Dynamic  TWTT  is  a  more  recent  development  that  allows  one  or  more  of  the  transceivers  to  be 
moving  [1,2]. 

Dynamic  TWTT 

Dynamic  TWTT  is  accomplished  in  the  same  fashion  as  static  TWTT,  with  the  exception  that  one  or  more 
of  the  receivers  is  moving.  The  moving  receiver(s)  introduce  motion-related  errors  that  are  not  present  in 
the  static  case.  A  dynamic  TWTT  configuration  is  illustrated  in  Figure  3. 

Not  all  of  the  cancellations  that  applied  to  the  static  TWTT  case  transfer  to  the  dynamic  case.  For  the 
dynamic  TWTT  scenario  in  Figure  3,  it  can  be  assumed  that  dAs  ~  dsA,  since  the  geostationary  satellite 
has  no  relative  motion  with  respect  to  the  earth  station  and  the  path  length  does  not  change.  This 
replicates  the  situation  in  the  static  case.  Unlike  the  static  case,  dSB  r  dus  for  the  dynamic  case,  since  the 
mobile  platform  has  moved  during  the  transmission  of  signals,  causing  the  transmit  and  receive  path 
lengths  to  be  different  between  the  geostationary  satellite  and  mobile  platform.  Because  the  mobile 
platform  is  in  motion,  the  Sagnac  effect  will  also  vary  and  produces  a  time-dependent  value. 


Figure  3.  Dynamic  TWTT  using  a  satellite  [2]. 


Taking  all  this  into  account,  the  time  differenced  measurement  for  dynamic  TWTT  becomes: 

A  -  B  =  —  [R(A)  -  R(B)  -  A PropDelay  +  SAb  -  SBa] 

where 


(1) 
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R(A)  =  time-interval  counter  reading  for  Station  A 
R(B)  =  time-interval  counter  reading  for  Station  B 
SAb  =  Sagnac  delay  from  Station  A  to  Station  B 
SBa  =  Sagnac  delay  from  Station  B  to  Station  A 
A  =  time  of  Clock  A 
B  =  time  of  Clock  B 

A PropDelay  =  change  in  propagation  delay  over  measurement  interval. 

The  A  PropDelay  term  is  a  time-varying  value  that  changes  based  on  the  relative  motion  of  the  mobile 
platform  as  well  as  how  the  velocity  vector  is  projected  onto  the  line  of  sight  vector  from  the 
geostationary  satellite.  The  Sagnac  delay  term  (Sab  -  Sba)  is  also  time-varying,  changing  based  on  the 
absolute  position  of  the  two  receivers  and  the  velocity  vector  projected  onto  the  equatorial  plane  [2], 


METHODOLOGY 

This  research  is  based  on  a  simulation  created  using  MATLAB*  and  contains  five  main  functions,  seen  in 
Figure  4.  The  load _params  function  involves  collecting  desired  input  parameters  from  the  user.  The 
generate  truth  function  uses  the  input  parameters  to  create  truth  data  that  will  simulate  the  environment 
that  is  being  measured.  The  generatemeas  function  uses  the  truth  data  to  generate  pseudorange  and 
TWTT  measurements  for  a  geostationary  satellite.  The  kalman  function  inputs  the  generated 
measurements  into  a  Kalman  filter  and  predicts  the  state  of  the  satellite  at  each  epoch  in  the  simulation. 
The  analyze  results  function  takes  the  results  of  the  Kalman  filter  and  compares  them  to  the  truth  data  to 
determine  the  accuracy  of  the  filter.  Each  of  these  functions  will  be  briefly  described  in  the  sections  that 
follow. 


'  ' 

load_params 


( - >1 

t  > 

r - \ 

generate_truth 

generate_meas 
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V _ J 

f - - 

analyze_results 

V _ > 


Figure  4.  Simulation  block  diagram 


Parameters 

The  simulation  begins  by  collecting  all  the  desired  input  values  for  a  host  of  variables  that  are  used 
throughout  the  simulation.  The  essential  parameters  are: 

•  Initial  ECI  state  for  the  geostationary  satellite 

•  Simulation  run  time  and  time  step  interval  (set  to  1  day  of  simulation  with  measurements  at  1- 
minute  intervals) 

•  GPS  satellite  ephemeris  date  selection 


514 


39th  Annual  Precise  Time  and  Time  Interval  (PTTI)  Meeting 


•  Clock  type  selection  and  clock  model  parameters 

•  Two-Way  Time  Transfer  measurement  accuracy 

•  Kalman  filter  parameters  (noise  values,  etc.) 

•  Monte  Carlo  parameters  (e.g.,  number  of  runs). 

Truth  Model 

The  truth  model  function  is  responsible  for  generating  all  data  that  will  be  considered  as  the  absolute 
truth.  There  are  three  main  operations  that  happen  in  this  function:  1)  propagation  of  the  geostationary 
satellite  state  forward  over  a  specified  time  interval,  2)  calculation  of  the  true  GPS  satellite  position  and 
clock  at  measurement  times,  and  3)  clock  error  modeling. 

GEO  Satellite  Propagation.  Implementation  of  a  simple  Kalman  filter  propagates  the  GEO  satellite 
state  vector  into  the  future.  The  initial  state  vector  x(ro)  is  provided,  along  with  an  initial  covariance 

matrix  P(ro)  and  a  dynamics  matrix  F(/0).  The  covariance  matrix  describes  the  accuracy  of  the  state 
vector  values  and  the  dynamics  matrix  explains  the  motion  of  the  state  vector. 

The  GEO  satellite  state  vector  is  comprised  of  three  position  and  three  velocity  states,  as  shown  in 
Equation  (2).  Note  that  the  simulation  is  implemented  within  an  Earth-centered  inertial  (EC-1)  frame,  and 
the  final  results  are  later  converted  to  an  Earth-centered  Earth-fixed  (ECEF)  frame  for  ease  of 
presentation  and  understanding. 


x  =  [X  Y  Z  X  Y  Z]T 

The  continuous  truth  model  for  propagating  the  GEO  orbit  is 

\(t)  =  F(/‘)x(t)  +  w(f) 

where  the  process  noise  w (t)  is  described  by 


(2) 

(3) 


£[w(f)w(f')]  =  Q(f  )<?(*-*') 


(4) 


The  dynamics  matrix  F(t)  describes  two-body  the  motion  of  the  GEO  satellite,  and  is  written  in 
continuous  form  as: 


F(0  = 


0 

0 

0 

_± 

r 

0 

0 


0  0  10  0 

0  0  0  1  0 

0  0  0  0  1 

0  0  0  0  0 

-4  o  ooo 

r 

0  -^-0  0  0 


(5) 


where  //  is  the  Earth’s  gravitational  parameter  and  r  is  the  current  orbital  radius  of  the  satellite. 
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A  real  satellite  orbit  would  include  the  effects  of  many  perturbing  forces,  such  as  higher-order  gravity 
terms,  the  effect  of  solar  pressure,  gravity  effects  from  the  moon  and  sun,  etc.  Many  of  these  effects  are 
well  understood  and  could  be  modeled  in  a  real  system.  These  types  of  deterministic  effects  were  not 
modeled  within  the  simulation,  because  they  would  basically  be  added  into  the  truth  model  and  then 
removed  in  the  filter  model,  having  no  significant  impact.  However,  it  is  understood  that  in  the  real 
world,  the  filter  model  will  not  perfectly  match  the  truth  model.  To  accurately  represent  this  effect  within 
the  simulation,  a  small  amount  of  dynamics  noise  was  added  to  the  truth  model  through  the  matrix  Q(t): 


Q(0  = 


0 

0 

0 

0 

0 

0 


0  0  0  0  0 

0  0  0  0  0 

0  0  0  0  0 

0  0  o-,,2  0  0 

0  0  0  a2  0 

0  0  0  0  a,,2 


(6) 


GPS  Satellite  True  Position  and  Clock  Calculation.  After  propagating  the  GEO  satellite  state  vector 
forward  in  time  over  the  entire  simulation  time  interval,  the  truth  model  function  then  calculates  the 
positions  and  clock  states  of  each  individual  GPS  satellite  over  the  entire  simulation  time  interval.  The 
true  GPS  satellite  position  and  clock  errors  were  calculated  using  a  precise  ephemeris  (.sp3)  file.  Later, 
within  the  estimation  Kalman  filter,  the  position  and  clock  error  are  calculated  using  the  broadcast 
ephemeris  data  for  the  same  time  period.  This  means  that  there  is  a  realistic  difference  between  the  truth 
and  the  modeled  satellite  positions  that  is  representative  of  a  particular  day  in  GPS  history. 


Clock  Modeling.  The  simulation  needed  to  model  the  receiver  clock  error  in  the  GPS  receiver  in  the 
GEO  satellite,  as  well  as  in  an  atomic  clock  on  the  ground  (used  in  one  of  the  scenarios  described  in  the 
results  section).  For  this  simulation,  the  clock  errors  were  modeled  using  a  clock  model  given  in  [7]. 


The  performance  of  atomic  clocks  was  simulated  using  a  three-state  polynomial  process  driven  by  white 
noise.  The  discrete  process  model  was  implemented  as: 
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where 


(7) 


xt  (tk )  and  xl  (tk+l )  =  clock  bias  error  at  times  tk  and  tk+l 
x,  (7, )  and  x2  ( tk  )  =  clock  drift  error  at  times  tk  and  tk+l 
x,  (7  )  and  x,  (7,  , )  =  clock  drift  rate  error  at  times  t,  and  t  , 
t  =  t,  -t=  time  interval 

/fc+1  k 

wj(k),  w3(k),  and  w3(k)  =  white  noises. 

The  clocks  cannot  be  modeled  deterministically  because  of  their  stochastic  nature.  Instead,  the 
performance  of  the  random  walk  noise  values  (vv/,  wj,  w3)  is  modeled  and  the  characteristic  Allan 
variance  curves  of  the  atomic  frequency  standards  are  matched  [6].  The  statistics  of  the  random  walk 
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noise  values  are  determined  by  the  values  of  the  variance  continuous  process  noise  power  spectral 
densities  (qi,  q2,  q3)  of  Qd  in  Equation  (8). 


Qd  (t)  =  E[w(k)w(k)T~\ 


qg  +  \q2P  +±qlT5 

2  8  ft1" 

iftr 

2  ft2" 

ftr  +  iftr3 

yftr: 

Us1 

1ft2"2 

ft2" 

(8) 


The  results  given  later  use  both  atomic  (rubidium)  clocks  and  ovenized  crystal  clocks  (representative  of  a 
good  crystal  clock  that  would  be  used  in  a  typical  spaceborne  GPS  receiver).  The  q  values  for  a 
spacebome  rubidium  clock  for  this  simulation  were  chosen  by  leveraging  research  conducted  in  the  Clock 
Improvement  Initiative  [7]  and  are  displayed  in  Table  1.  The  ovenized  crystal  clock  parameters  are 
obtained  from  [6].  To  calculate  a  clock’s  three-state  random  process  in  the  simulation,  initial  clock  bias, 
drift,  and  drift  rate  values  are  selected  from  Table  1  and  then  propagated  using  Equation  (7).  The  Qd  from 
Equation  (8)  was  used  to  generate  properly  correlated  wy,  wy,  and  w3  terms  using  a  U-D  factorization 
technique. 


Table  1.  Process  noise  values  for  simulated  clocks. 


Rubidium  Clock 

Ovenized 
Crystal  Clock 

qi  (bias) 

1.11  x  10'22  s2/s 

1.6  x  10"21  s2/s 

q2  (drift) 

2.22  x  10"32  s2/s3 

3.2  x  10"21  s2/s3 

q3  (drift 
rate) 

6.66  x  10"45  s2/s5 

0  s2/s5 

Generated  Measurements 

The  measurement  generation  function  is  responsible  for  creating  pseudorange  measurements  by  using  the 
information  supplied  by  the  truth  generation  function.  Pseudorange  values  are  normalized  range 
measurements  with  the  addition  of  errors  due  to  pseudorange  measurement  noise,  GPS  satellite  clock 
bias,  and  receiver  clock  bias.  The  pseudorange  equation  is: 

p  =  >/(*“'  -O2  +(/“' -y,.J+{zsa'-zrJ  +c5trec  -cStsa'  +oPR  (9) 

where 

xsat ,  ysat ,  zsal  =  true  ECEF  position  of  the  satellite  (m) 
xrec  ’  y'rec  ’  z,-ec  =  true  ECEF  position  of  the  receiver  (m) 

St  =  receiver  clock  bias  (s) 

8t~  =  satellite  clock  bias  (s) 
uPR  =  pseudorange  error  (m) 
c  =  speed  of  light  (m/s). 
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An  important  part  of  the  GPS  measurement  model  is  determining  exactly  which  GPS  satellites  are 
“visible”  to  the  GEO  GPS  receiver  at  any  point  in  time.  Visibility  is  determined  by  a  combination  of 
received  signal  strength  and  a  model  of  GPS  receiver  sensitivity. 

Received  Signal  Strength 

To  determine  the  strength  of  the  GPS  signal  that  is  received  by  the  GEO  satellite,  the  satellite  nadir  look 
angles  are  needed.  If  the  GEO  satellite  and  GPS  satellite  positions  are  known,  simple  vector  math  will 
produce  the  angles  6  and  a  (referenced  in  Figure  5),  which  are  the  GPS  satellite  look  angle  and  GEO 
satellite  look  angle,  respectively.  The  calculated  look  angles  are  then  used  with  antenna  gain  pattern 
information  to  determine  the  received  signal  strength.  Additionally,  any  signal  that  passed  within  400  km 
of  the  surface  of  the  Earth  was  deemed  unavailable  due  to  atmospheric  effects. 


Figure  5.  GPS  measurement  geometry. 


The  total  received  signal  power  can  be  calculated  according  to  [8]: 

P  _  PtGtGr  h2 

4 kR2  4 n 

where 


(10) 


Ps  =  received  signal  power  (watts) 

Pt  =  signal  power  at  transmit  antenna 
Gt  =  transmit  antenna  gain 
Gr  =  receive  antenna  gain 

R  =  distance  between  transmit  and  receive  antennas 
2  =  signal  wavelength  (GPS  LI  wavelength  «  5.255  m). 

The  resulting  value  is  the  signal  power  at  the  exit  of  the  receiver  antenna.  A  typical  transmit  (GPS 
satellite)  antenna  gain  pattern  was  obtained  from  [9].  The  GEO  satellite  receiver  was  modeled  using  the 
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gain  pattern  from  a  patch  antenna  that  flew  aboard  the  Falcon  Gold  experiment  from  the  Air  Force 
Academy  [10].  This  particular  antenna  is  representative  of  hardware  that  has  flown  on  previous  satellites. 

Once  the  total  receive  signal  power  was  determined,  a  C/N0  value  was  calculated  by  dividing  by  (or 
subtracting,  if  working  in  dB)  a  noise  power  density.  This  simulation  uses  a  standard  N0  value  of  202 
decibel-watts  [8]. 

GPS  Receiver  Model.  Once  the  received  signal  strength  has  been  calculated,  the  pseudorange 
measurement  noise  error  can  be  established  through  the  use  of  a  GPS  receiver  model.  Fundamentally,  the 
C/N0  value  defines  whether  or  not  a  pseudorange  measurement  is  available,  and  if  it  is,  what  the 
measurement  noise  would  be. 

Figure  6  shows  the  receiver  models  used  in  this  simulation.  The  best  performing  receiver  (labeled  “Fligh 
+”)  is  derived  from  data  obtained  in  [11],  and  it  is  essentially  an  optimal  receiver  which  performs  better 
than  most  real  receivers.  Similar  C/N0  data  is  located  in  [12].  For  each  simulated  receiver  model,  the 
point  at  which  the  line  ends  at  the  left  represents  the  minimum  C/N0  value  that  will  still  yield  a 
pseudorange  measurement. 


Two-Way  Time  Transfer  Measurements 

The  Two-Way  Time  Transfer  measurements  in  this  simulation  do  not  include  Sagnac  error  or  motion 
related  errors,  since  they  are  largely  deterministic  and  can  be  removed.  The  simulation  could  add  the 
errors  and  then  remove  them,  but  this  would  be  a  wasted  step  that  would  only  increase  computational  cost 
and  would  have  no  added  value.  For  simplicity,  this  simulation  assumes  that  the  propagation  delays  will 
cancel  as  in  the  static  TWTT  case.  The  resulting  TWTT  measurement  equation  is: 
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A  T=\  [TIC  (GEO)  -  TIC  (REF)]  +  uTWTT 

=  \\2cStCEO-lcStREF]  +  VTWTT 


C^^GEO  ^(itgEF  ^TWTT 


where 

StGEO  =  GEO  satellite  clock  error 
St REF  =  reference  clock  error 
v ^  =  TWTT  measurement  error 

TWTT 

c  =  speed  of  light. 


(11) 


TWTT  measurements  are  given  to  the  Kalman  filter  along  with  the  pseudoranges  for  measurement 
incorporation.  The  TWTT  measurement  noise  standard  deviation  values  used  in  this  simulation  are  10,  3, 
0.3,  0.03,  and  0.003  meters.  (Units  of  meters  rather  than  seconds  are  used  because  the  system  is 
ultimately  designed  to  navigate.  Dividing  these  values  by  the  speed  of  light  would  give  the  values  in 
seconds.) 


Kalman  Filter 

A  Kalman  filter  was  chosen  over  a  least-squares  batch  filter,  since  it  allows  for  the  use  of  new 
measurement  data  as  they  become  available  and  easily  models  stochastic  processes,  such  as  clock  errors. 
The  Kalman  filter  has  several  initial  values  that  govern  the  estimation  algorithm.  The  filter  must  know 
the  accuracy  level  of  the  incoming  measurements  and  how  much  to  trust  in  their  positioning  information, 
as  well  as  the  amount  of  process  noise  in  the  system. 

The  first  thing  needed  by  the  Kalman  filter  is  an  initial  state,  including  position,  velocity,  and  clock  error. 
These  initial  values  are  generated  by  adding  the  Gaussian  random  errors  shown  in  Table  2  to  the  truth 
data.  The  state  vector  is  defined  as: 


x  =  [x  Y  Z  X  Y  Z  cdtGE0  cSiGE0  cdtref  cSirefJ 

where 


(12) 


X,Y,Z  =  GEO  satellite  position  components 

X,Y,Z  =  GEO  satellite  velocity  components 

StGEO  and  SiGE0  =  GEO  satellite  clock  bias  and  clock  drift 

St  and  Si  ,  =  TWTT  reference  clock  bias  and  clock  drift 

ref  ref 

c  =  speed  of  light. 


Table  2.  Initial  state  error  standard  deviation  values. 


Initial  State 

Standard  Deviation  Value 

Position 

20  m 

Velocity 

0.01  m/s 

Clock  bias 

14  m 

Clock  drift 

20  m/s 
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Dynamics  Model.  The  dynamics  matrix  F  describes  the  motion  of  the  satellite,  and  is  generated  using 
the  same  equations  as  the  truth  generation.  The  Kalman  fdter  F  matrix  includes  the  clock  terms,  making 
it  a  10  x  10  matrix.  The  initial  covariance  matrix  P  describes  the  accuracy  of  the  state  vector,  and  will  be 
updated  as  the  filter  iterates.  The  process  noise  covariance  matrix  Q  describes  the  errors  associated  with 
propagating  the  state  covariance  matrix  P  through  time.  The  Q  matrix  includes  the  process  noise  value  of 
the  GEO  satellite  and  the  related  clock  q  values.  The  Q  matrix  is  equal  to  the  Q  matrix  used  in  the  truth 
generation  and  the  clock  q  values  are  taken  from  Table  1,  depending  upon  the  type  of  clocks  that  are  used. 
(The  only  difference  is  that  the  clock  and  position/velocity  estimation  are  all  done  simultaneously.)  The 
Q  matrix  used  in  the  Kalman  filter  does  not  change  throughout  the  simulation.  Standard  Kalman  filter 
propagation  equations  were  applied  within  the  simulation  [6]. 

Measurement  Model.  The  measurement  model  used  in  the  simulation  is 


z(t.)  =  h[x(t),t.]  + v(t) 

where 

z (t.)  =  measurement  vector  at  time  t,- 

\(t.)  =  zero-mean  white  Gaussian  measurement  noise  of  strength  R  =  E[vT(t,)  v(t,)] 
and,  for  GPS  measurements 


(13) 


h[xfr .),/,]  =  yjix" - O 2  +(/“' -yrJ2  +(z"  +C5tm  -cSts‘ 


The  measurement  equation  for  the  TWTT  measurements  was  given  in  Equation  (11). 


(14) 


After  the  h  vector  equations  are  written,  the  measurement  partial  derivative  matrix  El  is  constructed.  The 
H  matrix  relates  the  linearized  observations  to  the  estimated  states,  and  is  expressed  for  n  measurements 
as: 


where 


H 


H„ 


(15) 


H  _  Sh^ 

-,H2  = 

Sh2(x ) 

ft; 

ii 

Sx 

Sh,(x)  _ 

5h\{x) 

Sh,(x) 

Shx(x) 

Sx 

Sx j 

Sx2 

5X>n 

(16) 

(17) 


The  H  matrix  is  of  size  n  by  m,  where  n  is  the  number  of  measurements  and  m  is  the  number  of  states.  (In 
this  simulation,  m  =  10).  For  example,  if  there  are  two  pseudorange  measurements  and  one  TWTT 
measurement  at  a  given  epoch,  the  H  matrix  and  measurement  vector  z  will  be: 
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e'v 

< 

0 

0 

0 

1 

0 

0 

0 

H  = 

e; 

e; 

e2z  0  0  0 

1 

0 

0 

0 

0 

0 

0  0  0  0 

1 

0 

-1 

0 

A 

z  = 

Pi 

TWTT 

where 

e 

Xrec  ~  ^ 

e  —  rec 

X  //—  sat  ~  \2  /  sat \2  '  / 50/  \2 

V(X  ~xrec)  +(y  -yrec)  +  (Z  ~Zrec) 


y  - O2 + (/“'  - y,ec)2 + (*“  - O2 

sat 

Z  —Z 

e  =  rgc 

z  //  sat  \  2  .  /  \  2  .  /  sat  \  2 

V(X  ~xrec)  +(y  -yrec)  +  (Z  ~Zrec) 

In  the  example  given  above,  the  corresponding  R  matrix  is: 

0 

0 

2 

° TWTT  _ 

where 


R 


'PR  1 

0 

0 


0 

2 

7 PR  2 


(18) 


(19) 


(20) 


<JPR  =  standard  deviation  value  of  pseudorange  measurement  noise 
a twtt  =  standard  deviation  value  of  TWTT  measurement  noise. 

Standard  extended  Kalman  filter  measurement  incoiporation  equations  were  implemented  within  the 
kalman  function  of  the  simulation.  More  details  on  the  Kalman  filter  implementation  can  be  found  in 

[13]. 


RESULTS 

Analyzing  the  simulation  results  involves  comparing  the  Kalman  filter  estimated  state  with  the  true  state. 
The  analysis  depends  on  the  simulation  type,  being  either  a  single  run  or  a  Monte  Carlo  collection  of  runs. 
All  internal  calculations  within  the  simulator  were  performed  in  the  ECI  frame.  However,  prior  to 
interpreting  the  results,  the  errors  were  all  converted  to  ECEF  coordinates,  so  that  they  would  be  easier  to 
interpret.  This  is  mostly  due  to  the  fact  that  the  GEO  satellite  is  nearly  stationary  in  ECEF  coordinates. 

The  most  important  result  is  the  three-dimensional  positioning  error,  which  will  be  expressed  as  Mean 
Radial  Spherical  Error  (MRSE).  The  MRSE  is  analogous  to  a  three-dimensional  Distance  Root  Mean 
Square  (DRMS)  value.  For  a  Monte  Carlo  simulation,  the  MRSE  for  a  particular  epoch  is  calculated  by 
using  the  following  equation: 


522 


39th  Annual  Precise  Time  and  Time  Interval  (PTTI)  Meeting 


where 


MRSE  = 


£(.r2+ji2+z,2) 


(21) 


n  =  number  of  simulation  runs 

x.  =  X'",e  -  X[‘"er  =  difference  of  truth  and  fdter  X  position  for  epoch  i 
y i  =  Y'"‘e  -  T =  difference  of  truth  and  filter  Y  position  for  epoch  i 
z.  =  Z"”'  -  Zmer  =  difference  of  truth  and  filter  Z  position  for  epoch  i. 

The  MRSE  value  can  also  be  calculated  by  using  the  standard  deviation  values  that  exist  in  the  covariance 
matrix  for  a  particular  epoch,  as  seen  below. 


MRSE  = 


V  /  2  2  2  ' 

Lis  +y<  +z< , 


2  2 

(7  -her 


where 


(22) 


<jx  ,  cr  ,  a  =  standard  deviation  values  for  the  X,  Y,  Z  coordinates. 

Using  a  Monte  Carlo  simulation,  each  run  will  generate  different  position  values,  but  the  filter-computed 
covariance  values  will  be  the  same  for  every  single  run.  As  a  result,  the  covariance  standard  deviation 
values  from  a  single  run  can  replace  the  position  values  from  hundreds  of  runs  in  a  Monte  Carlo 
simulation,  if  the  filter  modeling  is  accurate. 

After  the  MRSE  is  calculated  for  each  time  epoch,  the  Root  Mean  Square  (RMS)  is  calculated  for  the 
entire  collection  of  epochs,  using  Equation  (23).  The  final  result  is  a  single  RMS  value  that  depicts  the 
level  of  error  in  the  estimated  GEO  satellite  position. 


where 


RMS  = 


n 


V  MRSE; 


(23) 


n  =  number  of  epochs 
MRSEj  =  MRSE  for  epoch  i. 

Covariance  Analysis  Validation 

As  a  first  step  in  validating  the  simulation  model,  several  Monte  Carlo  simulations  were  performed  to 
analyze  the  statistical  results  and  confirm  the  output  was  reasonable.  The  Monte  Carlo  simulations 
consisted  of  100  iterations.  Each  simulation  uses  the  same  parameters,  but  uses  a  different  set  of  random 
numbers  produced  by  the  random  number  generator  in  MATLAB  " .  Figure  7  shows  the  estimation  error 
in  the  Y  direction  for  100  Monte  Carlo  runs  using  a  “standard”  stand-alone  GPS  receiver  (i.e.,  no  TWTT 
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measurements).  The  areas  of  error  growth  occur  when  there  are  very  few  GPS  measurements  available. 
For  this  simulation,  the  GEO  was  placed  over  0°  latitude  and  0°  longitude,  so  the  Y  direction  represents 
the  local  vertical  direction. 

The  Monte  Carlo  ensemble  mean  and  standard  deviation  are  shown  in  Figure  8,  along  with  the  filter- 
computed  standard  deviation  for  the  same  runs.  Note  that  the  filter  does  a  good  job  of  capturing  the 
ensemble  results.  Other  scenarios  were  run  with  similar  correspondence  between  the  filter-computed 
values  and  the  Monte  Carlo  values.  Monte  Carlo  simulations  that  include  TWTT  measurements  had 
similar  performance  as  the  non-TWTT  case.  Based  upon  this  performance,  the  filter-computed  values 
were  deemed  to  be  sufficient  to  characterize  the  simulation  output.  By  doing  so,  results  were  able  to  be 
obtained  by  performing  a  single  run  of  the  filter  for  any  given  scenario. 
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Figure  7.  100  Monte  Carlo  runs,  standard  receiver,  ovenized  crystal  oscillator,  no  TWTT 
measurements. 


Scenario  Descriptions 

There  were  five  basic  scenarios  that  were  simulated: 

1.  Crystal  GEO,  no  TWTT.  This  scenario  is  considered  the  current  performance  baseline,  as  it 
represents  the  case  where  the  solution  is  based  entirely  on  a  single  GPS  receiver  in  orbit  on  the 
GEO  satellite  with  a  good-quality  ovenized  crystal  (non-atomic)  clock,  and  no  use  of  TWTT 
measurements. 

2.  Rb  GEO,  no  TWTT.  The  only  difference  between  this  scenario  and  Scenario  1  is  that  the  GPS 
receiver  in  orbit  is  now  driven  by  a  rubidium  clock  in  the  GEO  satellite.  There  is  still  no  TWTT 
involved. 

3.  Crystal  GEO,  Rb  TWTT.  In  this  case,  the  GEO  satellite  has  an  ovenized  crystal  clock,  but  a 
TWTT  system  is  used  to  measure  the  difference  between  the  onboard  clock  and  a  rubidium  clock 
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that  is  on  the  ground.  (TWTT  1  -o  accuracy  =  0.3m  ~  Ins).  This  TWTT  measurement  is  then 
incorporated  as  described  in  Equations  (18)-(20),  which  essentially  has  the  effect  of  correcting  the 
onboard  clock  to  match  the  ground-based  rubidium  clock. 

4.  Rb  GEO,  Rb  TWTT.  In  this  scenario,  there  is  both  a  rubidium  clock  onboard  the  GEO  satellite, 
and  a  TWTT  system  is  used  to  measure  the  difference  between  the  onboard  oscillator  and  a 
rubidium  clock  on  the  ground. 

5.  Crystal  GEO,  GPS  Time  TWTT.  In  this  final  scenario,  the  GEO  satellite  has  an  ovenized  crystal 
clock,  and  a  TWTT  system  is  used  to  measure  the  difference  between  the  onboard  clock  and  clock 
on  the  ground  that  is  slaved  to  GPS  system  time. 


Y-Direction  Position  Error  (ECEF) 


Figure  8.  Comparison  between  Monte  Carlo  ensemble  statistics  and  filter-computed 
statistics,  standard  receiver,  ovenized  crystal  oscillator,  no  TWTT  measurements. 


Simulation  Results 

Each  of  these  scenarios  was  run  for  each  of  the  six  levels  of  receiver  sensitivity  described  previously, 
using  a  TWTT  measurement  accuracy  of  Ins  (la).  Figure  9  shows  the  entire-run  RMS  values  for  each  of 
these  test  cases  (calculated  as  shown  in  Equation  23).  There  are  several  things  to  note  about  the  results 
shown  on  this  plot. 

First  of  all,  there  is  a  modest  performance  improvement  when  an  atomic  clock  is  used  rather  than  a  crystal 
oscillator.  Interestingly,  the  performance  was  almost  identical  among  Scenarios  2,  3,  and  4 — all  scenarios 
that  include  an  atomic  clock.  This  would  imply  that  using  a  TWTT  system  to  synchronize  to  an  atomic 
clock  on  the  ground  would  give  equivalent  navigation  performance  to  placing  an  atomic  clock  in  orbit. 
From  a  practical  point  of  view,  this  would  alleviate  the  need  to  place  atomic  clocks  in  orbit  for  many 
satellite  navigation  applications. 
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Secondly,  and  probably  most  importantly,  there  is  a  drastic  performance  improvement  when  the  TWTT 
system  is  used  to  synchronize  with  a  ground  clock  slaved  to  GPS  system  time  (Case  5).  The  overall 
improvement  in  positioning  performance  varies  between  60%-70%  across  all  receiver  models.  In 
essence,  this  approach  is  able  to  measure  the  absolute  error  in  the  GEO  GPS  receiver  clock.  As  a  result, 
each  GPS  pseudorange  measurement  actually  becomes  a  true  range  measurement.  (The  GPS  pseudorange 
measurements  are  called  pseudorange  measurements  rather  than  range  measurements  primarily  because 
they  include  the  effects  of  the  receiver  clock  error.)  Note  that  these  effects  all  happen  implicitly  within 
the  estimation  Kalman  filter,  as  it  estimates  both  the  receiver  clock  error  and  the  GEO  satellite  position 
simultaneously,  using  both  GPS  pseudorange  measurements  and  TWTT  measurements. 


GPS  Receiver  Sensitivity  Model 

Figure  9.  3D  RMS  position  error  vs.  GPS  receiver  sensitivity  levels  for  five  basic  scenarios. 


Finally,  it  is  important  to  point  out  that  the  relative  performance  improvements  (i.e.,  percentage 
improvement  over  Case  1)  are  similar  across  all  simulated  receiver  sensitivities.  This  implies  that  using 
TWTT  measurements  to  synchronize  with  a  GPS-disciplined  clock  is  of  value  in  many  different 
situations. 


It  is  insightful  to  evaluate  the  results  on  an  axis-by-axis  basis.  Figure  10  shows  the  RMS  1-a  errors  for 
both  Scenario  1  (no  TWTT)  and  Scenario  5  (TWTT  to  GPS-synchronized  ground  clock)  on  an  axis-by¬ 
axis  basis.  The  solid  lines  represent  Scenario  1,  and  the  dashed  lines  represent  Scenario  5.  It  is  clear 
from  this  plot  that  the  most  significant  benefit  of  using  the  TWTT  measurement  is  observed  in  the  X 
(vertical)  direction.  There  is  a  noticeable  improvement  in  the  Y  and  Z  directions,  particularly  for  the  less 
sensitive  receivers,  but  in  all  cases  it  is  the  X  direction  where  the  largest  improvement  can  be  seen.  This 
is  consistent  with  the  generally  understood  principle  from  GPS  that  the  receiver  clock  error  has  the 
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greatest  impact  in  the  vertical  direction.  It  makes  sense  that  using  a  TWTT  system  to  help  estimate 
receiver  clock  error  would  reduce  this  effect. 

Additional  Tests 

Variation  of  Ephemeris  Date.  To  ensure  the  simulation  does  not  depend  upon  the  ephemeris  data  for  a 
specific  day,  9  additional  days  were  selected  for  comparison.  One  day  was  selected  out  of  each  year  from 
1997  to  2006,  providing  a  comprehensive  evaluation  pool.  Each  simulation  was  identical,  other  than  the 
different  ephemeris  date,  and  used  the  worst-case  scenario  of  a  standard  sensitivity  receiver  with  no 
TWTT  measurements.  While  the  absolute  3-D  RMS  error  magnitudes  did  show  minor  variations  between 
different  ephemeris  sets  (up  to  17%  difference),  the  overall  trends  observed  in  Figure  10  were  nearly 
equivalent  for  all  ephemeris  sets.  It  is  apparent  that  differing  ephemeris  dates  do  not  have  a  significant 
impact  on  the  results  of  the  simulation. 
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Figure  10.  Axis-by-axis  evaluation  of  results  for  Scenarios  1  and  5. 


Required  accuracy  of  GPS  time  synchronization  for  Scenario  5.  In  Scenario  5,  the  TWTT 
measurement  is  taken  relative  to  a  clock  on  the  ground  that  is  slaved  to  GPS  system  time.  Inherently, 
there  will  be  errors  in  this  measurement  due  to  both  TWTT  measurement  accuracy  and  errors  in  the  GPS- 
slaved  clock  on  the  ground  (i.e.,  the  difference  between  the  reference  clock  and  true  GPS  time).  Both  of 
these  types  of  errors  will  have  the  same  effect  on  the  results,  so  they  are  considered  to  be  lumped  together 
as  a  total  TWTT  measurement  error  for  the  Scenario  5  cases. 

In  order  to  evaluate  the  sensitivity  of  the  results  to  this  total  TWTT  measurement  error,  the  Scenario  5 
results  were  recalculated  for  a  range  of  TWTT  measurement  error  values  between  0.3  ns  and  100  ns,  and 
the  results  are  shown  in  Figure  1 1  (along  with  the  baseline  case  with  no  TWTT  measurements,  for 
comparison  purposes).  It  is  interesting  to  note  that  the  simulated  results  are  not  highly  dependent  upon 
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TWTT  measurement  accuracy,  and  that  significant  degradations  in  performance  are  not  seen  until  the 
worst  case  value  of  100  ns. 


GPS  Receiver  Sensitivity  Model 

Figure  11.  3-D  RMS  position  error  for  varying  TWTT  measurement  errors  for  Scenario 
5  (TWTT  with  reference  clock  synchronized  with  GPS  Time). 


It  should  be  noted  that  all  errors  in  this  simulation  were  modeled  as  white,  Gaussian  errors,  and  that  if 
there  were  time-correlated  errors,  the  solutions  would  be  biased  and  not  be  as  good  as  the  simulation 
predicts.  (By  its  very  nature,  the  covariance  analysis  assumes  an  unbiased  solution).  Of  course,  the  level 
of  measurement  error  correlation  depends  upon  the  frequency  at  which  measurements  are  incorporated. 


CONCLUSIONS 

This  simulation  demonstrated  that  TWTT  measurements  can  be  extremely  useful  in  improving  the 
positioning  performance  of  high-altitude  satellites.  The  largest  benefit  comes  from  using  TWTT 
measurements  to  synchronize  between  the  GPS  receiver  clock  at  GEO  and  a  clock  on  the  ground  slaved 
to  GPS  system  time.  In  this  case,  3-D  RMS  positioning  accuracy  was  improved  between  60%-70%.  The 
level  measurement  error  of  the  TWTT  measurement  was  not  critical  to  the  results  for  this  simulation  as 
long  as  the  error  was  below  30  ns.  Also,  the  greatest  positioning  improvement  is  observed  in  the  vertical 
direction. 

A  more  modest,  but  still  significant  (2 1  %-38%),  performance  improvement  was  obtained  when  TWTT 
was  used  to  synchronize  between  the  GPS  receiver  at  GEO  and  an  atomic  clock  on  the  ground.  It  was 
shown  that  the  performance  in  this  case  was  equivalent  to  the  performance  of  placing  at  atomic  clock  in 
orbit. 
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Use  of  TWTT  measurements  allows  any  standard  GPS  receiver  to  operate  effectively  on  a  GEO  satellite 
with  reasonable  accuracy.  Accurate  GPS  navigation  in  high-altitude  orbits  provides  numerous 
opportunities,  such  as  automated  station-keeping  in  a  GEO  orbit.  Also,  by  substituting  automation  and 
removing  the  ground-based  ranging  systems,  the  cost  reduction  incurred  by  reducing  ground  support  is 
considerable. 


DISCLAIMER 

The  views  expressed  in  this  paper  are  those  of  the  authors  and  do  not  reflect  the  official  policy  or  position 
of  the  United  States  Air  Force,  Department  of  Defense,  or  the  U.S.  Government. 
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